function DecodingProbability_v6()

% clear workspace
clear all
close all
clc

% Configure parameters
T=100; % nb of transmitted packets
q=2; % field cardinality
k1=10; %
k2=10; %
k3=10; %
K1=k1; %
K2=k1+k2; %
K3=k1+k2+k3; %
g1=0.5; %
g2=0.5; %
g3=0.25; %
% n1= %
% n2= %
% n3= %

% Prepare decoding probability vectors
p_decode_l1=zeros(T,1);
p_decode_l2=zeros(T,1);
p_decode_l3=zeros(T,1);


% For layer 1


% For layer 2 - this might be right!!
for N=1:T % for number of transmissions
    for i=0:K1 % for each possible rank of layer 1 matrix
        val_1=ProbMatricesWithRank(g1*N,k1,i,q);
        val_2=ProbMatricesWithRank(g2*N,K2-i,K2-i,q);
        
        if g1*N<i
            val_1=0;
        end
        
        if g2*N<(K2-i)
            val_2=0;
        end
        N;
        i;
        val_1;
        val_2;
               
        p_decode_l2(N)=p_decode_l2(N)+val_1*val_2;
    end
end

p_decode_l2
length(p_decode_l2)

crap=0;
for w=1:T
    if p_decode_l2(w)==0
        crap=crap+1;
    else
        break;
    end
end
crap

% Plotting
figure(1)
hold('on')
% plot(1:T,p_decode_l1)
plot(1:T,p_decode_l2)
% plot(1:T,p_decode_l3)
hold('off')
grid('on')
pbaspect([2.5 1 1])
set(gca,'XTick',0:T/10:T)
xlim([0 T])
ylim([0 1])

% Save plot
print(gcf,'uep_ew_analytic.eps')

end


% Should work
function PMWR = ProbMatricesWithRank(m,n,r,q)

% Get first set of gaussian coefficients
gc=gausscoeffs(n,r,q);

% Calculate "sum"
val=0;
for k=0:r
    % This should be the one!
    val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(m*k+binomcoeffs(r-k,2)-n*m));
end

% Return probability of matrix 'm'x'n' with rank 'r'
PMWR=gc*val;

end

% Should work (Tested! see bottom)
function GC = gausscoeffs(m,r,q)
if r==0
    % disp('r = 0 in gauss coeffs')
    GC=1;
elseif r>0
    % disp('r > 0 in gauss coeffs')
    
    % Calculate numerator
    num=1;
    for w=m:-1:m-r+1
        num=num*(q^w-1);
    end
    
    % Calculate denominator
    denom=1;
    for w=r:-1:1
        denom=denom*(q^w-1);
    end
    
    % Calculate gaussian coefficient
    GC=num/denom;
    
elseif r<0
    disp('r < 0 error in gausscoeffs!!!')
end


end

% Not tested!!! ffs!
function bc = binomcoeffs(a,k)
% As on page 123 in "A course in combinatorics"

tmp_vector=ones(2,1);
tmp_index=1;

for w=0:-1:-k+1
    tmp_vector(tmp_index)=(a+w); 
    tmp_index=tmp_index+1;
end

num=prod(tmp_vector);

denom=factorial(k);
bc=num/denom;

end

function p = ksuccesntrials(n,k,p)
p=nchoosek(n,k)*p^k*(1-p)^(n-k);
end

















%% Testing Gaussian Coefficient generater

% Correct values are from: http://mathworld.wolfram.com/q-BinomialCoefficient.html
%
% m=2;
% r=1;
% q=2;
%
% if gausscoeffs(m,r,q) == (1+q)
%     disp('gauss test 1 succes!')
% end
%
% m=3;
% r=1;
% q=2;
%
% if gausscoeffs(m,r,q) == (1+q+q^2)
%     disp('gauss test 2 succes!')
% end






























% Deprecated (Numerical fix)
% Should work (No need to test)
% function DM = DifferentMatrices(m,n,q)
% DM=q^(m*n);
% end

%%  !!!!!!!!!!!!!!!!!junk!!!!!!!!!!!!!!!!111

%val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(n*k+nchoosek(r-k,2)));
% NOTE: We must not take nchoosek(0,2)!
% Are we making a mistake? The book says we should?
% Is this right?
% val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(n*k+nchoosek(r-k,2)));
